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Abstract 

We propose improvements in numerical evaluation of symmetric stable density and its par- 
tial derivatives with respect to the parameters. They are useful for more reliable evaluation of 
maximum likelihood estimator and its standard error. Numerical values of the Fisher information 
matrix of symmetric stable distributions are also given. Our improvements consist of modification 
of the method of Nolan (1997) for the boundary cases, i.e., in the tail and mode of the densities 
and in the neighborhood of the Cauchy and the normal distributions. 

1 Introduction 

There have been many researches to evaluate densities and quantiles of symmetric or general stable 
distributions. McCulloch (1998) considered efficient algorithms for approximating symmetric stable 
densities f{x;a) for the range a > 0.85, where parameter a denotes the characteristic exponent. 
Nolan (1997) gave accurate algorithms for general stable densities based on integral representations 
of the densities which were derived by Zolotarev (1986). Nolan provides a very useful program package 
"STABLE" on his web page*. However his program exhibits some unreliable behavior around the 
boundary as stated in the users guide of STABLE. Therefore even in the case of symmetric stable 
distributions, reliable computations of density functions including all the boundary cases is still 
needed. Furthermore for maximum likelihood estimation, it is desirable to directly compute the 
derivatives of the density function. In this paper we present reliable computations of symmetric 
stable density functions and their partial derivatives. Our computation of densities is accurate for all 
values of x and 0.1 < a < 2. Concerning the partial derivatives it is accurate in a somewhat smaller 
range of values. 

Regarding maximum likelihood estimation for the range a > 0.4, Nolan (2001) used interpolated 
stable densities and maximized the likelihood by approximate gradient search (constrained quasi- 
Newton method) because of its efficiency. But near the boundary of the parameter space interpolation 
may be inaccurate and the direct integral representation is used. Note that the direct integral 
representation is also not very reliable and slow near the boundary. In the symmetric case, Brorsen and 
Yang (1990) discussed maximum likelihood estimation using an integral representation of the densities 
given by Zolotarev (1986). But they have only considered the range a > 1 to avoid the discontinuity 
and nondifferentiability at a = 1. Furthermore they did not check the sample covariances of their 
maximum likelihood computation with the Fisher information matrix. 
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These previous researches on maximum likehhood estimation have not used direct evaluation of 
the derivatives of the log likelihood function with respect to the parameters. For reliable evaluation 
of the maximum likelihood estimator and its standard error, direct and reliable evaluation of the first 
and the second derivatives of the log likelihood function is desirable. 

The organization of this paper is as follows. In Section |21 we summarize notations and preliminary 
results on symmetric stable density. In Section |31 we provide an accurate algorithm for calculations 
of symmetric stable distributions which modifies Nolan (1997) for x near or oo and for a = 1 or 
a = 2 using various expansions. Accurate algorithms for the partial derivatives of symmetric stable 
distributions with respect to the parameters are given in Section and Section IHI Fisher information 
matrices are calculated in Section |H1 together with some simulation studies on the variance of the 
maximum likelihood estimator and the observed Fisher information. We also discuss behavior of 
Fisher information as q — > 2. Some discussions are given in Section [7| 



2 Preliminaries 

In this section we prepare notations and summarize preliminary results. There are many parameter- 
izations for stable distributions and much confusion has been caused. Our parameterizations follow 
the useful parameterizations for statistical inference which were given in Nolan (1998). 



2.1 Notations 

Let 

(2.1) ^{t) = ^{t;n,a,a) = ex.p{-\at\" + i^t) 

denote the characteristic function of symmetric stable distribution with parameters 

where a is the characteristic exponent, ^ is a location parameter and cr is a scale parameter. For the 
standard case {fJ-jCr) = (0, 1) we simply write the characteristic function as ^{t;a) = exp(— |t|"). 
The corresponding density is written as /(x; fi, a, a) and f{x; a) in the standard case: 

f{x;n,a,a) = -f{- — ^;a). 

a a 



At a = 1 



7r(l +x2) 
is the Cauchy density and at a = 2 

/(x;2) = -^exp(-xV4) 

is the normal density A^(0, 2). Accordingly the density can be defined to constitute a location-scale 
family. In the following, the first derivative of f{x;a) with respect to x is denoted by f'{x;a) and 
the second derivative is denoted by f"{x;a). Then 

(2.2) ^/(^; ") = «) = -^f'( ; 

ox on cr^ a 

(2.3) dx^^^^' ^' ^ 'djj?-^^^' ^' ^' ^ a^-^"^~a~'' 
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The partial derivatives with respect to a and a are written by subscripts, e.g., 

d d"^ 
fa{x; a, a) = Q^fix; cr, a), faa{x] fi, a, a) = g^fi^'^ «)• 

As above, when these derivatives are evaluated at the standard case {fi, a) = (0, 1) we write a), 
faa{x;a), etc. Note that 

(2.4) f^{x; a) = -f{x; a) - xf'{x; a), faa{x] a) = 2/(x; q) + 4x/'(x; a) + x'^f"{x; a). 
Furthermore we write ^ 

The reason we consider up to the second order derivatives of the density function is that in 
assessing the standard error of the maximum likelihood estimator 9, the observed Fisher information 

(2.5) W.,.....„)_--g_^^^ -nh\l^) '"h^(^ 

is usually preferred to the value of the Fisher information matrix at 9 (e.g. Efron and Hinkley (1978)). 

Note that there are other parameterizations of symmetric stable distributions than (|2.H) . However 
different parameterizations in the literature are smooth functions of each other including the boundary 
a = 2 and differentiations in terms of other parameterizations can be obtained from the results of 
this paper by the chain rule of differentiation. 

2.2 Preliminary results 

From equation (2.2.18) of Zolotarev (1986) or Theorem 1 of Nolan (1997), the density /(x; a) for the 
case a 7^ 1 and x > is written as 

7r 

(2.6) f{x;a) = — — / g{ip;a,x)ex.p{-g{ip]a,x))dip, 

TT\a — l\x Jq 

where 

, X cos 93^1 ~ cos(a - 1)(^ 

(2.7) g{v;a,x) - 
Note that at x = 



sm a(p J cos 93 



/(0;a) = -r(l + i 
for all < a < 2. 

For the case x — > and < a < 2, a 7^ 1, the following expansion can be used. 

Pc — k — 
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For a < 1, this series is not convergent but can be justified as an asymptotic expansion as x — > 0. 
For 1 < a < 2, it is convergent for every x. Similarly for the case x ^ oo and 0<a<2, we 
have 

J-/ \ 1 r(/ca + 1) . .-nak 

2.9 =-V^^n ^ -1 sinH^^ • 

k=l 

For a < 1 this series converges for every x ^ and for a > 1 this series can be justified as an 
asymptotic expansion as x — > oo. For a = 2 this asymptotic expansion is zero, which corresponds to 
the fact that the tail of normal distribution is exponentially small. These (asymptotic) expansions are 
stated in Bergstrom (1953), Section XVII. 6 of Feller (1971), Section 2.4 and Section 2.5 of Zolotarev 
(1986). 

3 Numerical evaluation of symmetric stable densities 

As in Nolan (1997) we numerically evaluate the density function using (|2.6j) . In (|2.6j) the func- 
tion g{ip; a, x) : [0, ^] — > [0, oo] plays an important role, because properties of this function make 
the numerical integration quite efficient. Note that g{ip; a, x) is continuous and positive, strictly 
increases from to oo for a < 1 and strictly decreases from oo to for a > 1. Therefore the inte- 
grand g{-) exp{—g{-)) is unimodal and its maximum value 1/e is uniquely attained at 931 satisfying 
g{(pi; a, x) = 1. When the value of density is small, the integrand concentrates around its mode very 
narrowly. Then quadrature algorithms may miss the integrand. Therefore we solve g{ip; a,x) = 1 for 
(fi and the integral is divided into two intervals around this mode (see Nolan (1997)). For numerical 
calculations of (|2.6() we use adaptive integration with singularity (QAGS) in GNU Scientific Library 
(2003). For most values of x and a this integration works well. 

However, when a is close to 2 this algorithm has some difficulty. Note that for a = 2 

g{v;2,x) = ( ^ ] 
\ 2 sm (/? y 

and g{Tr/2;2, x) = x/2. Therefore ipi exceed tt/2 when x > 2. There are some other numerical 
difficulties in 1)2. 6|) . We list these difficulties and propose alternative practical methods for evaluating 
the density. 

1. a is small and x ^ 0: 

If a is small, the density is very much concentrated at a; = 0. For example Nolan (1997) states 
/(0;0.1) = 1.155 X 10^ whereas /(0. 01; 0.1) = 1.66. In our calculations when a is small and x 
near the values of 1)2. 6() sometimes become larger than /(O; a), contradicting the unimodality 
of the stable density. For this case we can use the asymptotic expansion ()2.8|) . 

2. X — > 00: 

We cannot guarantee the accuracy of ()2.6|) in the case of x — > 00. Since stable distributions have 
heavy tails, reliable calculation of their densities is needed for large x. In our calculations when 
a > 1 and x is large the values of ()2.6|) sometimes become much smaller than the asymptotic 
expansion ()2.9() . For this case we can use the asymptotic expansion (|2.9|1 . 

3. a is near 1: 

The representation 1)2. 6|) can not be applied at q = 1 theoretically. The numerical quadrature 
of p. 6)1 becomes unreliable because of roundoff errors, when a is close to 1. 



4 



4. a is near 2: 

Though the representation (|2.6j) can be appUed near a = 2 theoretically, it seems to be too 
close to the normal distribution in the tail of the distribution. Actually the values of the density 
in the tail obtained by the integral representation (|2.6|) is much smaller than the asymptotic 
expansion (|2.9|) . 

For the rest of this section, we discuss the cases 3 and 4 above. 

For a = 1, we consider Taylor expansion of the density around a = 1. Let 

7 = 0.57722 

denote the Euler's constant throughout the rest of this paper. The Taylor expansion of /(x; a) around 
a = 1 is given as follows. 

(3.1) /(x; a) = fix; 1) + l)(a - 1) + l)(a - l)^ + l)(a - 1)=^ + o((a - l)^), 

2 D 

where 

1 ( x'^ — 1 / 1 \ 2x 

(3.2) /„(x;l) = (^l-7--log(l + x2)J +^^— -^arctanx 

(3.3) faa{x;i) = ~;;^YTx2)3~ I T V ~ ^ ~ 2 '-''^ ^ V -1-^'^'^*^"^ 

8x(x2-l) /3 1 2n 

H — 7— — arctan 2; - - 7 - ^ + 2; ) 



7r(l+x2)3""" \2 ' 2 

"*" 7r(l + x^)^ {'■''^ ~ "^^^"^ ~ ^ ~ ^ ~ ^'■''^ ^"^^ arctan x| . 

For convenience the explicit form of faaa{x; 1) is given in appendix 1)3.2(1 and (|3.3p are proved as 
follows. From the equation 4.40 on page 18 of Oberhettinger (1990) 

(3.4) / u'^-^e""" log u cos{uy)du = (a^ + y^y^Tii^) [cos(z^z){^(z^) - - log(o^ + y^)} - z sm{uz)] , 
Jo 2 

where 

Re > 0, z = arctan(y/a), = T' {i^)/T{h>). 

Differentiating (|3.4|) several times with respect to z^, setting 1/ = 2, and combining the results in the 
inversion formula we obtain (|3.2() and (|3.3p . The conditions for change of integral and differentiation 
are satisfied in these cases. Although higher order derivatives of /(x;q) with respect to a can be 
evaluated along the same line, we found that the three term expansion ((3.1)) is sufficiently accurate. 

Now we consider the case a — > 2 and x ^ 00. It seems natural to use asymptotic expansion 
(|2.9|) . However the normal density has an exponentially small tail and this expansion is meaningless 
for a = 2. However in view of smoothness at a = 2, an approximation around the normal density 
is desirable. This case is somewhat subtle, but we found that the following procedure works well 
numerically. Note that from (|2.8() . for each fixed x, f{x;a) is differentiable with respect to a (> 1) 
even at a = 2, i.e., for each fixed x we have 

(3.5) fix; a) = /(x; 2) + /,(x; a)(a - 2) + ^/^^(x; a)(a - 2f + o{{a - 2f). 
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Differentiating (|2.9|) witli respect to a, for large x, heuristically we have 

Ux;a){a - 2) ~ ^T{a + 1)^ cos (^) ^-"-^(a - 2) ~ ^(a + _ 

and 

/(x;a) ~ /„(x;a)(a - 2). 

We summarize our treatments of various boundary cases in Tabled The range of x and a and the 
number of terms k in the expansions H2.8() and (|2.9|) are shown. For 0.99 < a < 1.01 we use formula 
(|3.1|) and for a > 1.99999 and for x large we use the maximum of (|2.6|) and (|3.5j) . Note that for 
most cases a small number of terms in the expansions 1)2. 8() or ()2.9() is sufficient. The approximations 
are very effective since the expansions and the integral 1)2. 6p give virtually the same results for most 
values of x and a. 



Table 1: Approximations to stable density at boundary cases 



a\x 


X ^ 0, formula (ESI) 


X ^ oo, formula 1)2. 9|) 




number of terms k 


range of x 


number of terms k 


range of x 


[0.1, 0.2] 


k = 1 


X < 10"^^ 


k = W 


3 

X > 10 


(0.2, 0.5] 


k-5 


X < 10"^ 


(0.5, 0.99] 


X < 10"^ 


(0.99,1.01] 


formula (\'6.1\\ 


(1.01,1.99999] 


k = 10 


X < 10"^ 


k = 10 


X > 10i+« 


(1.99999,2.0] 


k = 85 


X <7 


max(((Tni),(E;SI)) 



4 Partial derivatives of the symmetric stable density with respect 
to the location and the scale parameters 

In this section we discuss the first and the second derivatives of the stable density with respect to the 
location parameter ji and the scale parameter a. 

4.1 The first derivatives w.r.t. location and scale 

Differentiating (|2.6() we obtain 

(4.1) f'{x;a) = --/(x;a) + " ^) / g{ip;a,x){l - g{ip;a,x))e^p{-g{ip;a,x))dip 

X Tlx y(y — 1) Jq 

for a ^ 1 and x > 0. In view of ()2.2|1 and 1)2.4(1 we only have to evaluate (|4.H) . The integrand 
g{-){l — g{-)) exp {—§{■)) has a positive local maximum and a negative local minimum. These are 
attained at ip2: g{x;a,<p2) = (local maximum) and (^3: g{x;a,(p3) = (local minimum). 

For a < 1 the ordering of these points is ^3 < fi < ^2 and the reverse holds for a > 1. Accordingly 
we divide the interval of integration and then the method of adaptive integration gives reliable values. 
Without these divisions the numerical integration sometimes does not converge or produces incorrect 
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values. Note that if a is very close to 2, the values of i = 1,2,3, are outside of the integration 
range. Therefore another treatment is needed for this case. 

As in the case of the density itself, we use alternative representations of the first derivative of the 
symmetric stable density around the boundary and following representations are needed. 

1. Expansions of the first derivative of the density. 



(4.2) =±fffl|±lM,_l)V-^ 



vra^ (2fc-l)! 

As in (|2.8jl . for 1 < a < 2, this series is convergent for every x. For a < 1, this series is justified 
as an asymptotic expansion as x ^ 0. Similarly for the case x — > oo we have 

(4.3) /'(x; a) = i f ^^^^(-1)'= sm(^)x---. 

^ k=\ 

For a < 1 this series converges for every x 7^ and for q > 1 this series is justified as an 
asymptotic expansion. 

2. Taylor expansion of derivative of density around a = 1: 

(4.4) /'(x; a) = f'{x; 1) + /^(x; l)(a - 1) + i/^ Jx; l)(a - 1)^ + o((a - 1)^), 



where 



fix; 1) = 2^ 



7r(l+x2)2^ 



and 



.// -.^ lf-2x3 + 6x/3 1, 2-6x2 



w . X 2x(x^ - 14x2 + 9) [7^2 / 1 \2 ' 

= V(l + x2)4 j Y+ (^1-^- 2^°g(^+^ -l-arctan^x 

8(3x^-8x2 + 1) /3 1 2^ 
arctan x I 7 log(l + x 



7r(l+x2)4 V2 2 
2x(x^ - 22x2 + 17) / 1 
vr(l + x2)4 (l-7-^log(l + x 

4(x'^ - 6x2 _^ g^(^2 _ 



arctan x + 



2U 



(l + x2)4 (1 + x' 

Table 1^ is a summary of approximations to the first derivative f'{x;a). The interpretation of 
Table[21is almost the same as Tabled Note that for 0.99 < a < 1.01 the approximation (|4.4j) around 
Cauchy using Taylor expansion to the order 0((a — 1)2) is accurate enough. In the range a S [0.1, 0.2), 
approximations are not good and we do not consider accurate calculations of derivatives. For a very 
close to 2 and x large, we have the same problem as in the case of the density. However ()4.1() gives 
reasonable values and we use (|4.1j) . Note that if a is very close to 2, then we do not observer very 
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Table 2: Approximations to f'{x;a) at boundary cases 



a\x 


X ^ 0, formula (|4.2|) 


X — > oo, formula l|4.3j) 




number of terms k 


range of x 


number of terms k 


range of x 


[0.2, 0.25] 


k = 5 


X < 10-« 


k = W 


3 

X > 10 


(0.25, 0.3] 


X < 10-6 


(0.3, 0.99] 


X < 10"^ 


(0.99,1.01] 


formula 


(1.01,1.99999] 


A; = 10 


X < 10-3 


k = 10 


x > 10i+'» 


(1.99999,2.0] 


formula (j4.ll) 




Figure 1: Derivative of density w.r.t. // 



Figure 2: Derivative of density w.r.t. a 



large x and the approximation to f'{x;a) is less important than the density itself. We will discuss 
this problem again in Section 6 in connection with accurate evaluation of the Fisher information for 
a very close to 2. 

We present the graphs of the derivative concerning // (Figure^, and the derivative concerning a 
(Figure 12) . Both derivatives are continuous in a and we do not see abrupt changes of the forms for 
different values of a. 



4.2 The second derivatives w.r.t. location and scale 

We only need to evaluate f"{x;a) as in the case of the first derivatives. Differentiating (|4.1|) we 
obtain 

(4.5) f"{x;a) = -^—-\f{x;a) + - — '^-f'{x;a) 

a — 1 x'^ a — I X 



a^sign(a — 1) I o, 

— , _ 3 / 9 {v;a,x){2 - g[(p-,a,x))eicp[-g[(p-,a,x))dip. 



Considering that f{x;a) and f'{x;a) can be calculated as in the previous sections, we only need to 
evaluate the integral which appear in the third term of (|4.5j) . The integrand g'^{-){2 — g{-)) exp(— g'(-)) 
in the second term on the right hand side of 1)4. 5|) has a zero at (^4: g{x; a, 994) = 2, a local minimum 
at ifi: g{x;a, ifi) = 1 and a local maximum at 995: g(x; 0,1/95) = 4. For a < 1 the order of these 
points is v^i < 934 < 935 and the reverse order holds for a > 1. Therefore we can divide the interval of 
integral according to these values and get accurate results. Alternative representations of the second 
derivative of the symmetric stable density around the boundary are as follows. 

1. Expansions of the second derivative of density. 



m. ^ {2k - 2)< 

For 1 < a < 2 this series is convergent for all x and for q < 1 this is justified as an asymptotic 
expansion as x — > 0. Similarly for x ^ 00 we have 

(4.7) fix; a) = 1 f; £M±^(_i)-i sm(^)x---3, 

k=l 

For a < 1 this series is convergent for x / and for 1 < a < 2 this series is justified as an 
asymptotic expansion. 

2. Taylor expansion of the second derivative of density around a = 1: 

(4.8) fix; a) = fix; 1) + fix; l)(a - 1) + ^/l(x; l)(a - 1)^ + o((a - 1)^), 



where 

fix; 1) 



Gx^ - 2 
7r(l + x2)3' 



w// ^1X^-6x2 + 1/11 ^ 3, 2x\ 12x(x2-l) 



9 



and 



= ^ ^(1^^2)5 ^ Y+ (1-7-2 log(l + x')j -l-^^^tan^x 



2(5x6-155x^ + 235x2-21)/^^ 1 , 2- 



WT^^ ['-'-2 + 

4x(llx^ - 70x2 _^ 3;^N 2(x6 - 50x^ + 85x2 

H 7 ;tte arctan x H ; tttb — 

7r(l+x2)5 7r(l + x2)5 



Table 3: Approximations to f"{x;a) at boundary cases 



a\x 


X ^ 0, formula (j4.6j) 


X ^ oo, formula (|4.7|) 




number of terms k 


range of x 


number of terms k 


range of x 


[0.2, 0.25] 


k = 5 


X < 10"^ 


k = 10 


3 

X > 10i+« 


(0.25, 0.3] 


X < lO"'' 


(0.3, 0.9] 


X < 10-^ 


(0.9, 0.99] 


X < 10"^ 


(0.99,1.01] 


formu 


a m 


(1.01,1.02] 


fc = 10 


X < 10~3 


(1.02,1.999] 


= 10 


3 

X > 10 1+" 


(1.999,2.0] 


formula (1131) 



Table 131 summarizes uses of various expansions. Here formula (|4.8|1 means use of Taylor expansion 
of the second derivative of density around a = 1 and 1)4. 5p means direct use of integral expression of 
the second derivative. Note for a ^ 2 and x ^ oo the same problem arise as in the density. However 
we do not use the second derivative of density for a £ (1.999, 2.0] in this paper. 

5 Partial derivatives with respect to the characteristic exponent 

In this section we discuss the first and the second derivatives of the stable density with respect to the 
characteristic exponent a. 

5.1 The first derivative w.r.t. a 

We only need to investigate the standard (/i = and a = 1) case. Although the representations of 
fa{x;a) sometimes become complicated, careful analysis of various cases gives accurate and efficient 
calculations. The direct differentiation of (|2.6j) yields the following representation. 

(5.1) /„(x;a) = — -/(x;a) + — — / ga{y:>;a,x){l - g{ip;a,x))exp{-g{ip;a,x))dip, 

a[l — a) TT\a — Ijx Jq 
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where 

(5.2) ga{x]a) 



xcosLp\ cos(a — 1) 



sin aip 

X 



cos If 



1 



(a - 1)2 
1 



log 



X COS if 



(a-1) 



■log 



sm aLf 
X cos 9? 



+ 



1 



a — 1 tan ac/? 



+ (/3tan(a — 1)^9 



sm aif 



+ 



= -9{x;a){hi{ip) + h2{ip) +hs{ip)). 
Here for convenience hi{(p), i = 1,2,3, are defined as 



1 



a — 1 tan aip 



+ (ptan{a — l)ip 



(5.3) 



hi{ip) 



1 



(a-1) 



■log 



X COS 9? 



sm a^j 



h2{ip) 



a — 1 tana(/9' 



/i3((^) = iptan{a — l)ip. 



We consider the integrand ga{-){l — g{-)) exp{—g{-)) in (|5.1j) separately, dividing it into three parts 

g{.)h,{-){l-g{-))eM-9{-)). 

For each g{-)hi{-){l—g{-)) exp{—g{-)), i = 1,2, 3, we divide the interval of integration appropriately 
based on zero points ipi, i = 1,4,5. The following properties are obtained by careful evaluation of 
the signs of derivatives of hi{-). 

hi{(p): strictly decreasing, zero at (p4: 



sin aip 



1, ifi < ip4 for a < 1 and ip^ < ipi for a > 1. 

h2{'-p)'- negative for a < 1. The sign changes from + to — at (/^s = ^ for a > 1. 

/i3((^): nonpositive for a < 1. nonnegative for a > 1. 

As in the density, we use alternative representations of the first derivative of the symmetric stable 
density concerning a around the boundary and following representations are needed. 

1. Expansions of the first derivative of the density w.r.t a. 



(5.4) 



{-x')K 



k=0 



{2k)\ 



For 1 < a < 2 this series is convergent for every x and for a < 1 it is justified as an asymptotic 
expansion as x — > 0. Similarly for x ^ oo we have 



/r r\ i- / \ l^r'(QA; + l), ^.k-\ ■ 

5.5 /„ x; a = "2^ /, ni ^™ 

I 1 T{ak + l) 
^vr^ {k-l)\ ^ ' 



vr 

— cos 
2 



~ka—l 



TTak 



logx sin 



nak 



-ka—l 



For a < 1 this series converges for every x ^ and for a > 1 this series is justified as an 
asymptotic expansion. 

2. Taylor expansion of the first derivative of a around a = 1: 

(5.6) fa{x; a) = fa{x; 1) + faa{x] l)(a - 1) + o(|q - 1|), 

where /q(x; 1) and faa{x; 1) are given by (|3.2j) and (|3.3|) . 
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In order to first confirm the derivative at a = 2 we utilize equation 3.21 on page 12 of Oberhettinger 
(1990). 



cos{uy)du 



1 

2 V2 



1 



(2a)- 



2 2 sec I —r- I exp 



'8a 



{D,iz)+D,{-z)). 



where z = {2a)~^y and Dy[z) is the parabolic cylinder functions (see p. 255 of Oberhettinger (1990)). 
Differentiating this with respect to setting u = 2 and adjusting some signs and constants we obtain 
fa{x\ 2). These values coincide with that of expansions in (|5.4j) at a = 2. Note the integral (|5.H1 is 
not accurate when a is very close to 2 especially for a > 1.9999999 . This inaccuracy of the integral 
(|5.1|) is illustrated in Figure El in comparison to the accurate values shown in Figure El 

Table 0] is a summary of approximations to the first derivative of the density function with respect 
to a. 

The formula 1)5.6(1 means Taylor expansion of the first derivative of a around a = 1. For a G 
(1.9999,2.0] we use only the expansion ()5.4j) and the asymptotic expansion ()5.5j) . 



Table 4: Approximations to fa{x]a) at boundary cases 



a\x 


X — > 0, formula EH) 


a; — > oo, formula(j5.5j) 




number of terms k 


range of x 


number of terms k 


range of x 


[0.1, 0.2] 


k = 1 


X < 10-^'^ 


k = 10 


3 

X > 10 i+" 


(0.2, 0.3] 


/fc = 5 


X < IQ-' 


(0.3, 0.5] 


X < 10"^ 


(0.5, 0.99] 


X < 10"^ 


(0.99,1.01] 


formu 


a (EH) 


(1.01,1.9999] 


= 10 


X < 10-^ 


k = 10 


3 

X > 10i+« 


(1.9999,2.0] 


= 85 


X < 8.0 


k = 20 


X > 8.0 



Here we present graphs of the first derivative w.r.t. a in Figure El Figure Figure El Figure El 
Figure 13 and Figure El The first derivative fa{x;a) with respect to a is shown in Figure El and El 
and the score function fa{x;a)/f{x;a) is shown in Figure |7| and El Note that graphs of fa{x;a) 
where a = 1.999999 and a = 2.0 in Figure El are almost the same corresponding to the fact fa{x; a) 
is continuous in a G (0, 2]. 



5.2 The second derivative w.r.t. a 

Differentiating ()5.H) we obtain 

2 



(5.7) faa{x;a) 
+ 



a(l — aY 
2sgn(a — 1) 



7r(a — lYx Jq 
a 



f{x;a) 

T 

' g{ip) {(1 + a)/ii(v?) + 2h2{ip) + hsiif)} (1 - g{ip)) exp{-g{ip))dip 



Tr\a — l\x 

a 

vrla — Ijx 



ai'^) [g (9?) - 35(V') + 1) {^i(v^) + ^2(V') + Hi^)} exp{-g{(p))dLp 







a Lp 
a — 1 sin^ aif cos2(a — l)ip 



(1 - 9{v>))ew{-9{v>))d^- 
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-0.05 

-0.1 - 
-0.15 - 

-0.2 
-0.25 

-0.3 
-0.35 




a=2.0 - 
a=1.5 - 
a=1.0 - 
a=0.8 




-0.15 
-0.2 

-0.25 
-0.3 

-0.35 



a=2.0 - 
a=1.5 - 
a=1.0 
a=0.8 

1.5 



Figure 3: First derivative fa{x]a) w.r.t. to a 



Figure 4: fa{x;a) (more detail around 0) 





Figure 5: fa{x]a) (accurate) 



Figure 6: fa{x;a) using (|5.1|) (not accurate) 




ct=2.0 
a=1.5 
a=1.0 
a=0.8 



1.5 



Figure 7: Score function fa{x]a)/f{x]a) Figure 8: fa{x]a)/f{x]a) (more detail around 0) 
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Concerning the second integral we utilize the arguments of Section 14.11 for efficient integration. The 
third integral is calculated by dividing the interval of integration using zeros of integrand (pi for 
i = 1,2,3. As for the fourth integral, we integrate separately the two terms in the integrand. An 
alternative representations of the second derivative of the symmetric stable density concerning a 
around the boundary are as follows. 

1. Expansions of the second derivative of density w.r.t. a. 



(5.8) 



2 ^ r((2fc + l)/a + l) 2^A: 

vtqS ^ {2k)\ ^ ' 



^ 00 

7ra^ 



r"{{2k + l)/a + 1) 



(2fc + l)(-x2)^ 



k=0 



{2k)\ 



For 1 < a < 2 this series is convergent for all x and for q < 1 this is justified as an asymptotic 
expansion as x — > 0. Similarly for x — > 00 we have 



(5.9) faa{x]a) 



1 ^T"{ak + l) 



vr 



(^-1) 



fc(-l 



sm 



Irak 



-ka—l 



k=l 

00 



k=l 



{k-l)l 

r(afc + 1) 
(^-1)! 



vr 



■ cos 



Irak 



k{-l 



,k~l 



log {x) - -r) sin 



logx sin 

Irak 



Irak 



-fco— 1 



vrlogx cos 



Tiak 



-ka—l 



For a < 1 this series is convergent for x 7^ and for 1 < a < 2 this series is justified as an 
asymptotic expansion. 



2. Taylor expansion of the second derivative of a around a = 1: 

(5.10) faa{x] a) = faaix; 1) + faaaix] l)(a - 1) + o{\a 

where faaa{x; 1) is given in Appendix 1X1 



Table El is a summary of approximations of the second derivative of density function concerning 
a. Formula (|5.1U|) means use of Taylor expansion of the second derivative of a around a = 1. For 
the sake of convenience we use the integral representation (|5.7j) when a £ (1.999,2.0) and x oo. 
This is somewhat problematic as in the case of the first derivative w.r.t. a and further investigation 
is needed for a near 2. In the calculations of the present paper, this problem appears only in the 
observed fisher information for a = 1.5. However the estimated values a very near a = 2 seldom 
occur and there seems to be no influence of this problem in the values presented in Table |H1 



6 Fisher informations of symmetric stable distributions 

Using the results of the previous section we can accurately evaluate the Fisher information matrices 
of the symmetric stable distributions. To the authors' knowledge there are only two somewhat 
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Table 5: Approximations to faaix',a) at boundary cases 



a\x 


a; — > 0, formula 1)5. 8(1 


X ^ oo, formula 1)5. 9|) 




number of terms k 


range of x 


number of terms k 


range of x 


[0.2, 0.3] 


k = 5 


X < 2.0 X 10- '' 


k = 10 


3 

X > 10 1+° 


(0.3, 0.5] 


X < 10-5 


(0.5, 0.99] 


X < 10--^ 


(0.99,1.03] 


formula 


(1.01,1.999] 


k = W 


X < 10-3 


k = 10 


3 

X > 10i+« 


(1.999,2.0) 


formula (|5.7j) 



incomplete results concerning the Fisher Informations of stable distributions. One is Nolan (2001) 
and the other is DuMouchel (1975). The former was based on numerical differentiations of densities 
and did not give the informations for a < 0.5. The latter had the problem of truncation of the integral 
and did not give the informations for a < 0.8. Considering that the tails and modes of distributions 
are important for calculating informations, our approach is desirable. We give the informations 
a G [0.2,2.0] of the symmetric distributions in TableEl In TablelHl'*' means that the quantity is not 
defined. The values in Table 101 largely coincide with those of DuMouchel (1975). 
The components of Fisher information matrix / are defined as follows. 

where (^i, 92,63) = (/i, cr, a). By symmetry it can be easily shown that I12 = I21 = 0. For notational 
clarity we write In = I^f_i, I22 = laa, hs = laa and 133 = laa- For numerical calculations of (|6.1j) we 
use adaptive integration of infinite intervals with singularities (QUAGIU) in GNU Scientific Library 
(2003). In can be shown that the information matrix at the Cauchy (a = 1.0) is analytically given 
as 



\2 



(6.2) I^^ = I,, = l/2a\ = ^(l-7-log2), = i | ^ + (7 + log 2 - 1) 

These values are consistent with numerical calculations in Table IHl 

In order to check our computations, we have done small simulation study of maximum likelihood 
estimation of a (for fixed fi = and o" = 1). We found that simulated information coincides well with 
the exact information. TableEJis the results of 1000 iterations of maximum likelihood estimation of a 
with the sample size of n = 50. We denote the maximum likelihood estimator by q. q is the average 
of 1000 iterations and is the variance of \/50 x (a — a). Though the convergence of 1/d"^ to laa 
seems somewhat slow for a < 1, the values of l/o"^ are consistent with the exact Fisher information 

laa ■ 

Furthermore we simulated the observed Fisher information ()2.5|) of a in Table |H1 Here Iaa{'^) 
corresponds to the right hand side of ()2.5() and 



1 / fa{xi;a)'^ ^ 
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Table 6: Fisher information matrix of symmetric stable distributions 



a \ lij 


T 




T 




T 


2.0 


0.5 




2.0 


OO 


* 


1.999 


0.4995 


1 


.9904 


29.461 


—0.8685 


1.99 


0.4960 


1 


.9321 


4.6197 


—0.6682 


1.95 


0.4842 


1 


.7631 


1.4108 


—0.4821 


1.9 


0.4727 


1 


.6127 


0.8846 


—0.3963 


1.8 


0.4552 


1 


.3898 


0.5937 


-0.3138 


1.7 


0.4424 


1 


.2189 


0.5028 


—0.2692 


1.6 


0.4334 


1 


.0775 


0.4726 


—0.2396 


1.5 


0.4281 





.9556 


0.4737 


—0.2174 


1.4 


0.4270 





.8475 


0.4973 


—0.1992 


1.3 


0.4310 


0, 


.7498 


0.5424 


—0.1832 


1.2 


r\ A A 1 r\ 

0.4419 





.6603 


0.6119 


—0.1679 


1.1 


0.4630 


0, 


.5774 


0.7132 


—0.1523 


1.05 


0.4790 





.5381 


0.7794 


—0.1440 


1.01 


0.4953 





.5075 


0.8413 


—0.1369 


1.0 


0.5 




0.5 


0.8590 


—0.1352 


0.99 


0.5049 





.4925 


0.8763 


—0.1332 


0.95 


0.5276 





.4631 


0.9552 


—0.1257 


0.9 


0.5641 





.4272 


1.0721 


—0.1154 


0.8 


0.6800 


0, 


.3586 


1.3928 


-0.0913 


0.7 


0.9094 


0, 


.2937 


1.8974 


-0.0611 


0.6 


1.4446 





.2325 


2.7414 


-0.0220 


0.5 


3.1167 





.1753 


4.2748 


0.0295 


0.4 


12.256 


0, 


.1226 


7.3994 


0.0979 


0.3 


188.09 


0, 


.0756 


14.925 


0.1869 


0.2 


149359.4 





.0367 


38.729 


0.2938 



Table 7: MLE (n = 50, 1000 iteration) 



a 


a 


1/O-Q 




1.99 


1.974 


4.4877 


4.6197 


1.8 


1.811 


0.6396 


0.5937 


1.7 


1.714 


0.5102 


0.5028 


1.5 


1.531 


0.4532 


0.4737 


1.3 


1.320 


0.5114 


0.5424 


1.0 


1.027 


0.7676 


0.8590 


0.8 


0.819 


1.1944 


1.3928 


0.5 


0.509 


3.8803 


4.2748 
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involves the first derivative only. Variance of /^^(l) and /aa(2) are also shown in the parentheses 
in Table IHl laaC^) has smaller variance than /^^(l) but some positive bias is observed in laaC^) 
compared to the true information laa- 



Table 8: Observed information (variance) (re = 50, 1000 iteration) 



a 


a 




(1) 




(2) 




1.5 
1.0 
0.5 


1.531 
1.025 
0.509 


0.4863 
0.8541 
4.2819 


(0.033) 
(0.157) 
(3.517) 


0.5145 
0.9001 
4.4660 


(0.018) 
(0.099) 
(2.755) 


0.4737 
0.8590 
4.2748 



For the rest of this section, we discuss behavior of the Fisher information as a t 2. Investigation 
of the Fisher information in this case is difficult because it involves detailed study of behavior of the 
score function as a t 2 and a; — > oo. 

DuMouchel (1975, 1983) have proved that 

laa — > cxD as a I 2. Nagaev and Shkol'nik (1988) made 

more detailed analysis of laa and shown that 

(6.3) ^- = 4Alog(l/A) (^ + --W)' ^ = '-"- 

Tableinigives the information for a extremely close to 2. /^^(l) means use of Taylor approximation 
around normal 1)3. 5|) and Icq, (2) means use of 1)2. 9() when x > 10 1+" for the density. The same is done 
with Jo-Q. There seems to be no substantial difference between /^^(l) and laaC^)- In the column N&S 
we show the values of l/(4Alog(l/A)). Our computation for A > 10~^ is accurate. The convergence 
of (|6.3|) seems to be very slow. Limiting theoretical behavior of I„a{ci) as a t 2 is not known at 
present. 



Table 9: Information around a = 2 



a 




Iaa{l) 


^aa (2) 


N&S 






2.0 - 10" 


-lu 


106860414 


92384764 


108573620 


-1.3482 


-1.3427 


2.0 - 10-^ 


-9 


10810787 


10167389 


12063736 


-1.3482 


-1.3395 


2.0 - 10' 


^8 


1144778 


1131645 


1357170 


-1.3482 


-1.3123 


2.0 - 10- 


-7 


127953 


127802 


155105.2 


-1.3478 


-1.2553 


2.0 - 10" 


-6 


14724 


14722 


18095.60 


-1.2094 


-1.1923 


1.99999 




1750 


1750 


2171.472 


-1.1100 


-1.1100 


1.9999 




217 


217 


271.4341 


-1.0069 


-1.0069 



7 Conclusions and some discussions 

In this paper we proposed reliable numerical calculations of the symmetric stable densities and their 
partial derivatives including various boundary cases. We found that except for very small values of a 
{a < 0.1) our method works very well. This enables us to reliably compute the maximum likelihood 
estimator of the symmetric stable distributions and its standard error. 
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For the family of stable distributions, the use of the observed Fisher information in H2.5() for 
assessing the standard deviation of the maximum likelihood estimator needs further investigation. 
Our simulation suggests that there may be some merit in using only the first term on the right hand 
side of (|23|) . 

Further study is needed to theoretically establish the limiting behavior of the Fisher information 
matrix as a t 2. 

Finally it is of interest to extend the methods of the present paper to the general asymmetric 
stable densities and to the multivariate symmetric stable densities. These extensions will be studied 
in our subsequent works. 



A The third order derivative of f{x;a) w.r.t. a around Cauchy 

Write 

poo 

G{u,y) = / ti'^'^e"" log^ ncos(ny)(in. 
Jo 

Differentiating (|3.4|1 three times and setting a = 1 we get 



■ log2(l + y2) _ r'(z.) log(l + y^) + r"(i.) 



1 



X ^ cos{uz) ( V'(^) ~ 2 ^°s(^ + ] ~ ^ sin(z^) 



+(1 + {-r{u) log(l + 2/2) + 2T'{u)} 

X I —z sin(i^z) ^'(/'(i^) — - log(l + y^)^ — z"^ cos(z^2;) + cos{i'z)ip' (vz) 

+{i + y^yl-T{u 



-z^ cos(z^z) ( V(i^) - 2 log(l + V^) 



-2z s\ii{y z)iIj' [u) + z^ smiyz) + cos{v z)il)" {u) 



where z = arctany. Then 



faaa{x; 1) = - (-G(4, x) + 3G(3, x) - G(2, x)) . 
vr 



B Derivatives of the density at a = 0.5 

For the purpose of checking some of our calculations, we can use the explicit formula of the density 
at the special case of a = 0.5. From (2.8.30) of Zolotarev (1986) /(a;;0.5) is written as 



/(x;0.5) 



_ 3 
X 2 



4x712 VV2^ 



+ cos,-LWl-cf^ 



where 



S{x) = sin {^t^^ and C{x) = 



COS ( -r 
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are known as Fresnel integral functions. We differentiate the above representation and obtain 



/'(x;0.5) 



and 

f"{x;0.5) 



3 x-2 

_ 7 
X 2 



4v/27r 



sm 



cos 



1 

4a; 
4x 



2 VV2^ 



2 



_ 7 
X 2 



4V27r 



15 



sm 



+ 



5x2 



4 V2^ 



cos 



4x 



4x 



5 



1 



^/2 

1 
2 

1 



5 



V2 



vrx 



+ ^°^'4^ 



sm 



4x 



V2 



vrx 



sm 



4x 



+ cos 



1 

4x 
C 



1 



\/27rx 
1 



V2 

1 

2 

1 



TTX 



c 



V2 



vrx 



+ 



1 



47r 



V27fx 
9x-^ 



vr 



We have confirmed that our formulas in Section numerically coincide with these explicit expressions 
at a = 1/2 including the boundary cases. 
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